Phase transition of a two dimensional binary spreading model 
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We investigated the phase transition behavior of a binary spreading process in two dimensions for 
different particle diffusion strengths (D). We found that N > 2 cluster mean-field approximations 
must be considered to get consistent singular behavior. The TV = 3, 4 approximations result in 
a continuous phase transition belonging to a single universality class along the D € (0, f ) phase 
transition line. Large scale simulations of the particle density confirmed mean-field scaling behavior 
with logarithmic corrections. This is interpreted as numerical evidence supporting that the upper 
critical dimension in this model is d c — 2. The pair density scales in a similar way but with an 
additional logarithmic factor to the order parameter. At the D = endpoint of the transition line 
we found DP criticality. 



I. INTRODUCTION 

The study of nonequilibrium phase transitions in sys- 
tems with absorbing phases is an active area of research 
in Statistical Physics with applications in various other 
fields such as chemistry, biology and social sciences [||. 
The classification of the types of phase transitions found 
in these systems into universality classes is, nevertheless, 
still an incomplete task. 

The directed percolation (DP) universality class is the 
most common nonequilibrium universality class . Di- 
rected percolation was indeed found to describe the criti- 
cal behavior of a wide range of systems, despite the differ- 
ences in their microscopic dynamic rules. However, the 
presence of some conservation laws and/or symmetries in 
the dynamics has been found to lead to other universality 
classes Q. 

The pair contact process (PCP) || is one of the models 
whose (steady state) critical properties belong to the DP 
universality class. If the model is generalized to include 
single particle diffusion (PCPD or annihilation/fission 
model ) , a qualitatively distinct situation arises since 
states with only isolated particles are no longer frozen 
and the question was raised of whether this would mod- 
ify the universality class. A field theoretical study of 
the annihilation/fission model was presented long ago 0. 
Unfortunately, it relies on a perturbative renormalization 
group analysis which breaks down in spatial dimensions 
d < 2 so that the active phase and the phase transition 
are inaccessible to this study. The upper critical dimen- 
sion d c is 2 for this bosonic theory, where multiple site 
occupancy is allowed, contrary to the usual lattice models 
and Monte Carlo simulations. A fermionic field theory 
is not available but it is expected to lead to d c = 1 g - 
therefore mean field predictions, with some logarithmic 
corrections, would be seen in d = 1 if the latter is the 
correct theory. 

Monte Carlo, Coherent Anomaly [p|JlQ[ and DMRG 
studies §Q of the 1 - d PCPD proved to be rather 



hard due to very long relaxation times and important 
corrections to scaling. Several hypothesis were put for- 
ward in order to classify its critical behavior: single type 
flfTll versus two regions of different behavior parity 
conserving (PC) class ||, mean field behavior, diffusion- 
dependent exponents. Some related models were also 
studied (lf-[l|] with the aim of identifying the relevant 
features which determine the critical properties. The 
matter is not yet fully clarified, but it seems more likely 
that this system belongs to a distinct, not previously en- 
countered, universality class. Park et al |l5|] have also 
pointed out that the binary character of the particle cre- 
ation mechanism, rather than parity conservation, might 
be the crucial factor determining the type of critical be- 
havior of PCPD. Higher dimensional studies of PCPD- 
like models are thus necessary in order to clarify their 
universal properties and thus contribute to a full under- 
standing of the nonequilibrium phase transitions puzzle. 

In the present work we have studied a two-dimensional 
model where particle diffusion competes with binary cre- 
ation and annihilation of pairs of particles. The model 
is described in the following section. Cluster mean field 
studies are presented in section III and Monte Carlo sim- 
ulations are discussed in section IV. Finally we summa- 
rize and discuss our results. 



II. THE MODEL 

The sites of a square lattice of side L are either oc- 
cupied by a particle (1) or empty (0). The following 
dynamic rules are then performed sequentially. A par- 
ticle is selected at random and i) with probability D is 
moved to a (randomly chosen) empty neighbor site; with 
the complementary probability, and provided the parti- 
cle has at least one occupied neighbor, then ii) the two 
particles annihilate with probability p or iii) with proba- 
bility 1—p two particles are added at vacant neighbors of 
the initial particle. The selection of neighbors is always 
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done with equal probabilities; the updating is aborted 
and another particle is selected, if the chosen sites are not 
empty/occupied as required by the process. In reaction- 
diffusion language one has 
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It is clear that, in the absence of particle diffusion (D = 
0), only sites which belong to a pair of occupied neighbors 
are active. In terms of pairs - taken as redefined 'par- 
ticles' - we have a unary process where 'particles' are 
destroyed at rate p or give birth to new 'particles' at rate 
1— p; the number of offsprings is greater or equal to two - 
because new pairs may also be formed with next nearest 
neighbor particles - and parity of the number of 'parti- 
cles' is not conserved. This is similar to the PCP and 
one expects to see a phase transition, in the DP univer- 
sality class, between an active phase with a finite density 
of pairs (at low p) and an absorbing phase without pairs 
(for p> p c ). 

When particle diffusion is included, one has a qualita- 
tively different situation, since configurations with only 
lonely particles are no longer absorbing - the only ab- 
sorbing states are the empty lattice and the configura- 
tions with a single particle. There is parity conservation 
in terms of particles and the creation and annihilation 
mechanisms are binary. The nature of the phase transi- 
tion, expected to occur at some value p c (D), is investi- 
gated below. 



III. CLUSTER MEAN-FIELD CALCULATIONS 



We performed TV-cluster mean-field calculations [ |l 9|j2 0| 
for this model. Since the details of the dynamics will not 
influence the values of the mean field critical indices, we 
have considered a simpler one dimensional version of the 
model where the creation takes place at the nearest and 
next-nearest sites to one side of the mother particle. 

At the site (N = 1) level, the evolution of the particle 
density p (denoted by n in &) can be expressed as 



dp 

at 



-2pp 2 + 2(1 -p)p 2 ! (1 -pf 
which has a stationary solution 
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with p c — 1/2. The pair density P2 (c in Q notation) 
is just the square of p at this level. For p <~ p c the 
densities behave as 



p(oo) cx [p c - p) 13 
p 2 (oo) oc {pc-p) P2 



(4) 
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with (3 = 1 and (32 
the critical point 



2 leading order singularities. At 
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which implies the leading order scaling is 

p(t)<xt- a , p 2 (t)Kt- a * , (7) 
with a = 1/2 and a 2 = 1, while in the absorbing phase 

ptfioct- 1 , p 2 (t)(xt- 2 . (8) 

All these exponents coincide with those found for the 
PCPD model || at the same level of approximation. 

In the pair (N = 2) approximation, the density of 'I' p 
and the TP pair density p 2 are independent quantities. 
One can easily check that the evolution of particles can 
be expressed as 

dp 

at 

l-2p + p 2 



= -2p(l - D)p 2 + 
+ 2(l-D)(l-p)p 2 (p 
while the evolution of pairs 
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+ (l-D)(l-p)p 2 (p-p 2 )(l-2p + p 2 ) i>j2 

Owing to the nonlinearities, we could not solve these 
equations analytically and had to look for numerical so- 
lutions. The critical indices thus obtained at different 
diffusion rates D are shown in Table ffl. As we can see, 
there are two distinct regions. For D >^ 0.2 p c is con- 
stant and (3-2 = 2, while for D <~ 0.2 p c varies with D 
and (3 2 = 1. All these results are in complete agreement 
with those of the PCPD model in the pair approximation. 

In the N = 3 level approximation the situation 
changes, as we can see in Table |: the two distinct re- 
gions for D > disappear and (3 2 = 2 everywhere as 
found in the site approximation. At D = 0, however, 
the particle density does not vanish at the transition but 
goes to p(p c ) = 0.2931. This means that the N = 3 
level approximation is already capable to describe the 
absorbing state that contains frozen, isolated particles. 
For p <~ p c , p(p c ) — p oc (p c — p)@ with (3=1, the same 
critical exponent as the order parameter (the pair den- 
sity) therefore we redefine eq.(4) now. These results are 
also in agreement with those of the PCP model |2 
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TABLE I. Summary of N — 2, 3, 4 approximation results 



2 



This kind of singular mean-field behavior persists for 
N = 4 (Table | Fig.0) and can be found in the N = 3, 4 
level approximations of the PCPD model as well p3[ . 
These results suggest that the N = 2 approximation is 
an odd one. Recent investigations in similar PCP-likc 
models |24 25 have also shown discrepancies in the sin- 
gular behavior of the low-level cluster mean-field approx- 
imations. One can conclude that in these models at least 
N > 2 cluster approximations are necessary to find a cor- 
rect mean- field behavior. It is probably just a coincidence 
that the N = 1 calculation produced the same results. 
We shall thereafter ignore the N = 2 results and refer to 
the N — 1,3 scenario as the mean field prediction. 
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FIG. 1. N = 4 cluster mean-field results for p2- The curves 
correspond to D = 0,0.01,0.05,0.1,0.2,0.3,0.4,0.5,0.75 (left to 
right). The same kind of convex curvature corresponding to 
fa = 2 can be observed for D > 0, while it is different for 
D = (dashed line) corresponding to (3 = 1. 



IV. SIMULATION RESULTS 

The simulations were started on small lattice sizes 
(L = 100, 200) to locate the phase transition point 
roughly at D = 0,0.05,0.2,0.5,0.8. The particle den- 
sity decay p L (t) was measured up to t max = 60000 MCS 
in systems started from fully occupied lattices and pos- 
sessing periodic boundary conditions. Throughout the 
whole paper t is measured in units of Monte-Carlo sweeps 
(MCS). For D = 0.05 we have not done so detailed anal- 
ysis as for other diffusion rates but only checked that 
the results are in agreement with the conclusions derived 
from the D = 0.2, 0.5, 0.8 data. 

Then we continued our survey on larger lattices: L = 
400, 500, 1000, 2000 and determined p c at each size by 
analyzing the local slopes of p(t) 

-\n[p{t)/p(t/m)] 



eff(t) = 



(11) 



ln(m) 

(we used m = 8). In the t -> oo limit the critical curve 
goes to exponent a by a straight line, while sub(super)- 
critical curves veer down(up) respectively. The p c (L) 



estimates exhibit an increase with L hence at the true 
critical point the critical like a e f t curves of a given L are 
sub-critical. This excludes the possibility of finite size 
scaling study at p c . 



A. Dynamical scaling for D > 

For the largest system size (L — 2000) at D = 0.5 
diffusion rate the local slopes analysis results in a p c = 
0.43915(1) and the corresponding a = 0.50(2) decay ex- 
ponent agrees with the mean-field value (see Fig.@). 
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FIG. 2. Local slopes of the particle density decay at 
D — 0.5 and L = 2000. Different curves correspond to 
p = 0.4392, 0.43916, 0.43913, 0.4391 (from bottom to top). 

We also measured the pair density p2 (t) ; applying a lo- 
cal slope analysis similar to eq.(pd|) suggests (Fig.||) the 
lack of phase transition of this quantity at p = 0.43915. 
Instead the curves veer up which may lead to different 
p c and CV2 estimates. Such strange behavior has already 
been observed in PCPD model simulations 711281. 
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FIG. 3. The same as Fig.|| for p 2 (t). 
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An explanation for this discrepancy was pointed out by 
Grassberger in the case of the PCPD model Ran- 
dom walks in two dimensions are just barely recurrent 
and single particles can diffuse very long before they en- 
counter other particles. Therefore it is natural to expect 
that p(t)/p 2 (t) ~ hi(i) and Fig.^j shows this really hap- 
pens at p c . 
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FIG. 4. p(t)/p2(t) and logarithmic fit for critical curves 
determined by p(t) analysis. The top curves correspond to 
D = 0.8, p = 0.475, the middle ones to D = 0.2, p = 0.41235, 
while the bottom ones correspond to D — 0, 5, p = 0.43916. 
Note that the ratio is smallest at D = 0.5. 

Therefore at D — 0.5 we can conclude that a 2 ~ a ~ 
0.5 taking into account logarithmic corrections. This 
however contradicts the mean-field approximation value 
a 2 = 1. 

Similar local slopes analysis for D = 0.2 and D = 0.8 
seem to imply a = 0.46(2) and a = 0.57(2) respectively. 
First this raises the idea that the exponents would change 
continuously with D as it was observed in some one di- 
mensional PCPD simulations |p|Jl8|]. Nevertheless the de- 
viations from 0.5 are small hence we tried to fit our data 
including logarithmic corrections. Logarithmic correc- 
tions may really arise if d c = 2 as predicted by bosonic 
field theory 0. The precise form of these corrections is 
however not known for the present case so we have tried 
several functional dependences and found that 

((a + bln(t))/t) a (12) 

is a good choice. As Fig.|] shows for D — 0.2 this really 
works with a — 0.507 exponent. 

Similarly for D = 0.8 the same logarithmic formula 
fitting resulted in a = 0.497. The coefficient of the loga- 
rithmic correction term is negative (b = —0.2776), while 
it is positive for D = 0.2 (b = 0.468). 
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FIG. 5. Logarithmic fit (circles) for p(t) (upper curve) and 
with the form f jl3| ) (diamonds) for /92(f) (lower curve) at 
D = 0.2 and p = 0.41235. 

These results suggest that logarithmic corrections to 
scaling should work for all cases we investigated but at 
D = 0.5 they are very small and change sign. Indeed ap- 
plying the same formula for the D = 0.5, p = 0.43913 
data we obtained a = 0.496 with b = 0.00027 and 
a = 1.552. As we found logarithmic corrections to the 
particle density decay and a logarithmic relation between 
P2(t) and p(t) we may expect even stronger logarithmic 
corrections to the p^{t) data. Trying different forms for 
D — 0.2, 0.5, 0.8 we found that taking into account m 2 (f) 
correction terms is really necessary and the best choice 
is 

((a + b\n{t) + cln 2 (t))t)- a . (13) 

This resulted in a 2 = 0.5007 for D = 0.2 (see Fig]|); 
a = 0.501 for D = 0.5 and a = 0.484 for D = 0.8. 
All these results imply that a — a 2 = 0.5 independently 
from the diffusion rate D. For p(t) this agrees with the 
mean-field approximations and we don't see a change of 
universality by varying D inferred from the N — 2 ap- 
proximation. The critical behavior of p 2 {t) however dif- 
fers from the a 2 IF — 1 prediction. 

B. Static behavior for D > 

The p c estimates for different sizes were used to ex- 
trapolate to the true critical value. Simple linear fitting 
as a function of l/L resulted in the values given in Ta- 
ble ||. For determining steady state exponents the densi- 
ties p L (t,p, D) and p\(t,p, D) were followed in the active 
phase until level-off values were found to be stable. Av- 
eraging was done in the level-off region for 100 — 1000 
surviving samples - those with more than one particle 
jl2| . Again at each p and D we extrapolated as a func- 
tion of l/L to the linif^oo p L (oo,p, D) values. The local 
slope analysis of exponent (3 
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a (c ln[p(oo,e. t ,£i)] - ln[p(oo, e^-i, £>)] , . 

p e ff{ei,V) = — — , (14) 

ln(c<) - ln(ej_i) 



(where e = p c — p) shows that the order parameter p 
exhibits a (3 ~ 1 asymptotic scaling at D = 0.2,0.5,0.8 
(Fig.[|) although a correction to scaling can be seen in all 
cases. 
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FIG. 6. Effective order parameter exponent results. Cir- 
cles correspond to D = 0.2, diamonds to D = 0.5, squares to 
D = 0.8 data. 
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FIG. 7. Logarithmic fitting to £2(00) at D = 0.5 using the 
form eq.(p^|). The coefficients are a = 0.112, b = 0.01 and 
P = 0.96(5). 



C. Data collapse for D > 

To test further the possibility of mean-field critical be- 
havior we performed finite size scaling on our p L (t,p, D) 
data assuming the mean-field exponents ]7| /? = 1, v j_ = 1 
and the scaling form 



(16) 



As Fig.g shows, a good data collapse was obtained for 
p c = 0.4395(1) at D = 0.5. 



The /3 ~ 1 agrees with the mean-field prediction. Do- 
ing the same analysis for p2 the local slopes seem to ex- 
trapolate to (3 ~ 1.2 for each £>. This is very far from the 
mean-field value j3*f F = 2 and we don't see any change 
by varying the diffusion rate down to D = 0.05. We 
have investigated the possibility of different logarithmic 
corrections and found that the 



p=(e/(a + bHe)r (15) 



form gives very good fitting with j3 = 0.96(5) for D = 0.5 
while for D = 0.2 and D = 0.8 (similarly to the exponent 
Q!2 case) we need to take into account hi 2 (e) correction 
terms to obtain similar good fitting (see Table || and 
Fig.^). Therefore we concluded that, as in a case, the 
steady state exponents are equal: j3 = /3a- Note that we 
have checked that logarithmic corrections to scaling can 
also be detected in p data with f3 ~ 1 exponent. 
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FIG. 8. Finite size data collapse according to scaling form 
© for L = 200,400,500,1000,2000. Different symbols denote 
data for p = 0.436,0.437,0.4375,0.438,0.4382,0.4384,0.4386 

Similarly the scaling form 

p L (t,e,D)<xt-P/^g(te^) (17) 



5 



(a = f3/v\\) can be checked near p c , assuming the mean- 
field values [0/3=1 and v\\ = 2. For the largest size 
(L = 2000) at D = 0.5 the best collapse of curves cor- 
responding top= 0.438,0.4382,0.4384,0.4386,0.4388 was 
obtained for p c — 0.4394(1). This agrees well with previ- 
ous p c estimates within margin of numerical accuracy. 
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sity decay. This process was solved exactly by Lee |2!| 
who predicted the following late time scaling behavior in 
d= 2 



Pit) 



1 



8irD 



\nt/t + 0{l/t) 



(19) 



We measured p(t) at p = 0.45 and D = 0.5 in a L = 2000 
system up to t max = 3 x 10 5 MCS. As Figure [To| shows, 
for intermediate times the density decays faster than this 
power-law in agreement with results for PCPD |28| but 
later crosses over to the expected eq.(|l9|) behavior with 
amplitude 0.078(2) and with a 4.46/f correction to scal- 
ing term. 

Unlike what we found at the critical point (see sec- 
tion B), p2(t) decays faster than pit) in the absorbing 
phase. The long-time behavior seems to be p2 oc t~ 2 
which agrees with the mean-field prediction. 
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FIG. 9. Data collapse according to scaling form ( |17| ) 
at D = 0.5. Different curves correspond to data at 
p = 0.438,0.4382,0.4384,0.4386,0.4388 



D. The D = case 

As explained above, we expect that this model exhibits 
2+1 dimensional DP universality because for the pair 
density the conditions of the DP hypothesis (||| are sat- 
isfied. Indeed at p c = 0.3709(1) we found that the de- 
cay exponent of pairs is 0.1 = 0.45(1) and the steady 
state density approaches zero with the scaling exponent 
/?2 = 0.582(1) in agreement with the estimates for this 
class a = 0.4505(10) and (3 = 0.583(14) At the crit- 
ical point the density of isolated particles takes a nonzero 
value, usually called the natural density, p{p c ) — 0.135. 
In p(| we showed that in case of the PCP and an other 
Id model exhibiting infinitely many absorbing states the 
nonorder field follows the scaling of the order parame- 
ter field. Here we found that the total density shows a 
singular behavior 

pip)- p(p c )<x{p c -pf (18) 

with the redefined exponent = 0.60(2) agreeing with 
that of the DP class within margin of numerical accu- 
racy. 
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FIG. 10. Density decay in the inactive phase. For large 
times we found p(t) = (4.46+0.078(2) ]nt)/t behavior (dashed 
line) by fitting data. 
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TABLE II. Summary of simulation results at criticality 



E. Scaling in the inactive phase 

According to the bosonic field theory [Q in the inactive 
phase the A + A — » reaction governs the particle den- 
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V. CONCLUSIONS 



We have investigated the phase transition of a two 
dimensional binary spreading model exhibiting parity 
conservation. In what concerns cluster mean-field ap- 
proaches, the results are similar to those of the PCPD 
model at the corresponding level of approximation |6|,p3[ . 
The N — 2 results suggest two different universality 
classes depending on the diffusion strength. Higher 
(N = 3,4) order cluster mean-field show a single univer- 
sality class characterized by f3 = 1 and a — 1/2. Com- 
paring these with other recent results for PCP-like mod- 
els and with the simulations, we believe that the N = 2 
case yields spurious results — although two universality 
classes were apparently observed in a study of the one 
dimensional PCPD || — and so N > 2 cluster approx- 
imations are necessary to describe the mean-field singu- 
larity correctly. This is not surprising and was already 
found in similar models p4] , p5[ . Note that in both the 
N = 3 and N = 4 approximations the p c seems to have 
a discontinuity by approaching D = 0. Similar disconti- 
nuity in the phase space of the PCP model was recently 
reported as the result of an external particle source J25| . 
This behavior may be the subject of further studies. 

We performed extensive and detailed simulations along 
the phase transition line and found a single universal- 
ity class with the order parameter exponents (3 = 1 and 
a = 0.5 for all D > 0. Logarithmic corrections to scaling 
were detected that are weakest at D = 0.5. In the lack 
of a theoretical prediction, we have selected the best log- 
arithmic fitting forms taking into account up to 0(ln 2 ) 
terms, but we cannot rule out the possibility of other 
logarithmic correction forms. Scaling function analysis 



confirmed the i/± — 2 and v\ | — 1 mean- field values. This 
seems to indicate that the critical dimension is d c = 2 
as predicted by the bosonic field theory. In the inactive 
region, the decay of particle density at larg e times was 
found to agree with an exact prediction |B9| 

The pair density p2 for p < p c (where the bosonic 
field theory breaks down) was shown to exhibit the same 
singular behavior as the order parameter, apart from a 
logarithmic ratio. Simulation results of the PCPD model 
p8| , fTo[ found indications for similar behavior. The rea- 
son why the mean-field approximation fails to describe 
the singular behavior of p2 is yet not clear to us but in 
the two-component description of the model it indicates 
strong coupling between pairs and particles (similarly to 
other models [p|]3l|]). In the inactive region, however, p2 
and p scale differently. At the D = endpoint of the 
transition line we found 2+1 dimensional DP critical be- 
havior of p2 with infinitely many frozen absorbing states 
similarly to the PCP model. 

We have found identical predictions for the present and 
the PCPD model within mean-field, which seem to be 
confirmed by our simulations - and also by preliminary 
simulations for the 2 - d PCPD (27||2§]. One thus con- 
cludes that it is very likely that parity conservation is 



irrelevant for this transition, as in the one-dimensional 
case |H| and in certain models with exclusion [§2| . Fur- 
ther renormalization group studies of these systems are 
necessary for a proper justification of these results. 
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